<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of demspeech_dual</title>
  <meta name="keywords" content="demspeech_dual">
  <meta name="description" content="DEMSPEECH_DUAL  Sigma-Point Kalman Filter based Speech Enhancement Demonstration.">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../../menu.html">Home</a> &gt;  <a href="#">ReBEL-0.2.7</a> &gt; <a href="#">examples</a> &gt; <a href="#">dual_estimation</a> &gt; demspeech_dual.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../../menu.html"><img alt="<" border="0" src="../../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="menu.html">Index for .\ReBEL-0.2.7\examples\dual_estimation&nbsp;<img alt=">" border="0" src="../../../right.png"></a></td></tr></table>-->

<h1>demspeech_dual
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>DEMSPEECH_DUAL  Sigma-Point Kalman Filter based Speech Enhancement Demonstration.</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>This is a script file. </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre class="comment">  DEMSPEECH_DUAL  Sigma-Point Kalman Filter based Speech Enhancement Demonstration.

  A single phoneme of speech, corrupted by additive colored noise is enhanced
  (cleaned up) through Dual SPKF (SRCDKF) based estimation.

  A single speech phoneme sampled at 8kHz is corrupted by additive colored (pink)
  noise. We use a simple linear autoregressive model (10th order) to model the
  generative model of the speech signal. We model the pink noise by a known 6th
  order linear autoregressive process driven by white Gaussian noise with known
  variance. The SNR of the noisy signal (y=clean+noise) is 0dB.

  The colored noise modeling (augmented state space model) is done according to
  the method proposed in: &quot;Filtering of Colored Noise for Speech Enhancment and
  Coding&quot;, by J. D. Gibson, B. Koo and S. D. Gray, IEEE Transactions on Signal
  Processing, Vol. 39, No. 8, August 1991.

  See also : GSSM_SPEECH_LINEAR
   Copyright (c) Oregon Health &amp; Science University (2006)

   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for
   academic use only (see included license file) and can be obtained from
   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the
   software should contact rebel@csee.ogi.edu for commercial licensing information.

   See LICENSE (which should be part of the main toolkit distribution) for more
   detail.</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../../matlabicon.gif)">
<li><a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>	ADDRELPATH  Add a relative path which gets expanded into a absolute path</li><li><a href="../../.././ReBEL-0.2.7/core/geninfds.html" class="code" title="function InferenceDS = geninfds(ArgDS)">geninfds</a>	GENINFDS  Generate inference data structure from a generalized state space model and user defined inference parameters.</li><li><a href="../../.././ReBEL-0.2.7/core/gennoiseds.html" class="code" title="function NoiseDS = gennoiseds(ArgDS)">gennoiseds</a>	GENNOISEDS    Generates a NoiseDS data structure describing a noise source.</li><li><a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>	GENSYSNOISEDS  Generate process and observation noise data structures for a given InferenceDS data structure</li><li><a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>	REMRELPATH  Remove a relative path (which gets expanded into a absolute path)</li><li><a href="../../.././ReBEL-0.2.7/core/srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>	SRCDKF  Square Root Central Difference Kalman Filter (Sigma-Point Kalman Filter variant)</li><li><a href="demspeech_dual.html" class="code" title="">demspeech_dual</a>	DEMSPEECH_DUAL  Sigma-Point Kalman Filter based Speech Enhancement Demonstration.</li><li><a href="../../.././ReBEL-0.2.7/examples/gssm/gssm_speech_linear.html" class="code" title="function [varargout] = model_interface(func, varargin)">gssm_speech_linear</a>	GSSM_SPEECH  Generalized state space model for single phoneme speech enhancement</li></ul>
This function is called by:
<ul style="list-style-image:url(../../../matlabicon.gif)">
<li><a href="demspeech_dual.html" class="code" title="">demspeech_dual</a>	DEMSPEECH_DUAL  Sigma-Point Kalman Filter based Speech Enhancement Demonstration.</li></ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre>0001 <span class="comment">%  DEMSPEECH_DUAL  Sigma-Point Kalman Filter based Speech Enhancement Demonstration.</span>
0002 <span class="comment">%</span>
0003 <span class="comment">%  A single phoneme of speech, corrupted by additive colored noise is enhanced</span>
0004 <span class="comment">%  (cleaned up) through Dual SPKF (SRCDKF) based estimation.</span>
0005 <span class="comment">%</span>
0006 <span class="comment">%  A single speech phoneme sampled at 8kHz is corrupted by additive colored (pink)</span>
0007 <span class="comment">%  noise. We use a simple linear autoregressive model (10th order) to model the</span>
0008 <span class="comment">%  generative model of the speech signal. We model the pink noise by a known 6th</span>
0009 <span class="comment">%  order linear autoregressive process driven by white Gaussian noise with known</span>
0010 <span class="comment">%  variance. The SNR of the noisy signal (y=clean+noise) is 0dB.</span>
0011 <span class="comment">%</span>
0012 <span class="comment">%  The colored noise modeling (augmented state space model) is done according to</span>
0013 <span class="comment">%  the method proposed in: &quot;Filtering of Colored Noise for Speech Enhancment and</span>
0014 <span class="comment">%  Coding&quot;, by J. D. Gibson, B. Koo and S. D. Gray, IEEE Transactions on Signal</span>
0015 <span class="comment">%  Processing, Vol. 39, No. 8, August 1991.</span>
0016 <span class="comment">%</span>
0017 <span class="comment">%  See also : GSSM_SPEECH_LINEAR</span>
0018 <span class="comment">%   Copyright (c) Oregon Health &amp; Science University (2006)</span>
0019 <span class="comment">%</span>
0020 <span class="comment">%   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for</span>
0021 <span class="comment">%   academic use only (see included license file) and can be obtained from</span>
0022 <span class="comment">%   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the</span>
0023 <span class="comment">%   software should contact rebel@csee.ogi.edu for commercial licensing information.</span>
0024 <span class="comment">%</span>
0025 <span class="comment">%   See LICENSE (which should be part of the main toolkit distribution) for more</span>
0026 <span class="comment">%   detail.</span>
0027 
0028 
0029 <span class="comment">%===============================================================================================</span>
0030 
0031 clear all; close all; clc;
0032 
0033 <span class="keyword">if</span> ~exist(<span class="string">'aryule'</span>)
0034   error(<span class="string">' [demspeech_dual] This demonstration requires the Matlab Signal Processing Toolbox to function correctly.'</span>);
0035 <span class="keyword">end</span>
0036 
0037 help <a href="demspeech_dual.html" class="code" title="">demspeech_dual</a>
0038 
0039 disp(<span class="string">' '</span>);
0040 disp(<span class="string">' '</span>);
0041 
0042 disp(<span class="string">'Two speech time-series estimates are extracted from the estimated state vectors.'</span>);
0043 disp(<span class="string">'The first is generated by taking the first component (zero''th lag term) of the'</span>);
0044 disp(<span class="string">'state vector. The second estimate is generated by using the last (10th lag)'</span>);
0045 disp(<span class="string">'component of the state vector, which is a fixed-lag smoothed estimate (it uses'</span>);
0046 disp(<span class="string">'more data).'</span>);
0047 disp(<span class="string">' '</span>);
0048 disp(<span class="string">'After each iteration (over the whole speech sequence) of the filter, the normalised'</span>);
0049 disp(<span class="string">'MSE of each estimate is displayed. Three speech sequences are also played over the'</span>);
0050 disp(<span class="string">'audio device: The first is the noisy sequence, the second is the first estimate and'</span>);
0051 disp(<span class="string">'the third is the second (full-lag) estimate.'</span>)
0052 disp(<span class="string">' '</span>);
0053 
0054 dosound = input(<span class="string">'Do you want to enable the audio component of this demo (0=no 1=yes) ? '</span>);
0055 
0056 <span class="comment">%--- General setup</span>
0057 
0058 <a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>(<span class="string">'../gssm'</span>);         <span class="comment">% add relative search path to example GSSM files to MATLABPATH</span>
0059 <a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>(<span class="string">'../data'</span>);         <span class="comment">% add relative search path to example data files to MATLABPATH</span>
0060 
0061 
0062 <span class="comment">%-- Load clean speech, noise and noisy speech (0dB SNR)</span>
0063 
0064 load speech_data;                       <span class="comment">%</span>
0065 
0066 B=0;
0067 N=1500;
0068 
0069 clean = clean(B+1:B+N);
0070 noisy = noisy(B+1:B+N);
0071 noise = noise(B+1:B+N);
0072 
0073 <span class="comment">%-- Display speech waveforms</span>
0074 
0075 figure(1);clf; subplot(<span class="string">'position'</span>,[0.025 0.1 0.95 0.8]);
0076 p1=plot(noisy,<span class="string">'k+'</span>); hold on;
0077 p2=plot(clean,<span class="string">'b'</span>);
0078 xlabel(<span class="string">'time'</span>);
0079 legend([p1 p2],<span class="string">'noisy'</span>,<span class="string">'clean'</span>);
0080 title(<span class="string">'ReBEL Speech Enhancement Demo'</span>);
0081 axis tight
0082 drawnow
0083 
0084 <span class="comment">%-- Initialise GSSM data structure</span>
0085 
0086 model = <a href="../../.././ReBEL-0.2.7/examples/gssm/gssm_speech_linear.html" class="code" title="function [varargout] = model_interface(func, varargin)">gssm_speech_linear</a>(<span class="string">'init'</span>);            <span class="comment">% initialize</span>
0087 
0088 <span class="comment">%=====================================================================</span>
0089 <span class="comment">%  Generate InferenceDS data structures for dual estiamtion. We need</span>
0090 <span class="comment">%  one for the state estimator and one for the parameter estimator.</span>
0091 
0092 ftype = <span class="string">'srcdkf'</span>;                                     <span class="comment">% we will use square-root central difference Kalman filter (SRCDKF)</span>
0093                                                       <span class="comment">% estimator</span>
0094 
0095 paramParamIdxVec = 1:model.speech_taps;               <span class="comment">% index vector of the system parameters to be estimated (don't estimate</span>
0096                                                       <span class="comment">% colored noise model parameters)</span>
0097 
0098   <span class="comment">%-- Setup state estimator</span>
0099   Arg.type = <span class="string">'state'</span>;                                   <span class="comment">% inference type (state estimation)</span>
0100   Arg.tag = <span class="string">'State estimation for GSSM_SPEECH system.'</span>; <span class="comment">% arbitrary ID tag</span>
0101   Arg.model = model;                                    <span class="comment">% GSSM data structure of external system</span>
0102 
0103   InfDS_SE = <a href="../../.././ReBEL-0.2.7/core/geninfds.html" class="code" title="function InferenceDS = geninfds(ArgDS)">geninfds</a>(Arg);                             <span class="comment">% Create inference data structure and</span>
0104   [pNoise_SE, oNoise_SE, InfDS_SE] = <a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>(InfDS_SE, ftype);   <span class="comment">% generate process and observation</span>
0105                                                              <span class="comment">% noise sources for state estimator</span>
0106 
0107   <span class="comment">%-- Setup parameter estimator</span>
0108   clear Arg;
0109   Arg.type = <span class="string">'parameter'</span>;                                <span class="comment">% inference type (parameter estimation)</span>
0110   Arg.tag = <span class="string">'Parameter estimation for GSSM_SPEECH system.'</span>; <span class="comment">% arbitrary ID tag</span>
0111   Arg.paramFunSelect = <span class="string">'both'</span>;                           <span class="comment">% We use the full system dynamics as observation, i.e. obs=hfun(ffun(x))</span>
0112   Arg.paramParamIdxVec = paramParamIdxVec;               <span class="comment">% parameters to be estimated index vector (don't estimate colored</span>
0113                                                          <span class="comment">% noise model)</span>
0114   Arg.model = model;                                     <span class="comment">% GSSM data structure of external system</span>
0115 
0116   <span class="comment">%-- Explicitely define a observation noise source for the parameter estimator. This is needed for the colored noise</span>
0117   <span class="comment">%   case, since it uses an implicit (within the augmented state) observation noise formulation. When a parameter estimator</span>
0118   <span class="comment">%   is derived from this type of model, one has to override the default (empty/dummy) observation noise source that is</span>
0119   <span class="comment">%   generated.</span>
0120 
0121   Arg.model.Ndim = 1;                                    <span class="comment">% We need to override these field for the parameter estimator</span>
0122 
0123   oNoise_Arg.type = <span class="string">'gaussian'</span>;                          <span class="comment">% and actually define a true observation noise source.</span>
0124   oNoise_Arg.cov_type = <span class="string">'full'</span>;
0125   oNoise_Arg.dim = 1;
0126   oNoise_Arg.mu = 0;
0127   oNoise_Arg.cov  = sqrt(pNoise_SE.cov(2,2));
0128 
0129   Arg.model.oNoise = <a href="../../.././ReBEL-0.2.7/core/gennoiseds.html" class="code" title="function NoiseDS = gennoiseds(ArgDS)">gennoiseds</a>(oNoise_Arg);
0130 
0131   InfDS_PE = <a href="../../.././ReBEL-0.2.7/core/geninfds.html" class="code" title="function InferenceDS = geninfds(ArgDS)">geninfds</a>(Arg);
0132 
0133   [pNoise_PE, oNoise_PE, InfDS_PE] = <a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>(InfDS_PE, ftype);    <span class="comment">% generate process and observation</span>
0134                                                               <span class="comment">% noise sources for state estimator</span>
0135 
0136 
0137 <span class="comment">%-- ESTIMATE SIGNAL</span>
0138 
0139 N = length(noisy);                                    <span class="comment">% number of samples in frame</span>
0140 
0141 Xh_SE = zeros(InfDS_SE.statedim, N);                  <span class="comment">% setup estimation buffers</span>
0142 Xh_PE = zeros(InfDS_PE.statedim, N);                  <span class="comment">%     &quot;              &quot;</span>
0143 
0144 init_mod = aryule(noisy, model.speech_taps);          <span class="comment">% initial model is fit to noisy speech</span>
0145 init_mod = -1*init_mod(2:end);
0146 
0147 Xh_PE(:,1) = init_mod(:);                             <span class="comment">% initial model</span>
0148 
0149 Sx_SE = eye(InfDS_SE.statedim);                       <span class="comment">% initial Cholesky factor of SE estimate covariance</span>
0150 Sx_PE = eye(InfDS_PE.statedim);                       <span class="comment">% initial Cholesky factor of PE estimate covariance</span>
0151 
0152 InfDS_SE.spkfParams = [sqrt(3)];                      <span class="comment">% CDKF scale parameter for SE estimator</span>
0153 InfDS_PE.spkfParams = [sqrt(3)];                      <span class="comment">% CDKF scale parameter for PE estimator</span>
0154 
0155 number_of_runs = 10;                                  <span class="comment">% number of iterations over data</span>
0156 
0157 mse = zeros(2,number_of_runs);                        <span class="comment">% mean square error buffer</span>
0158 
0159 pNoise_PE.cov = 1*eye(InfDS_PE.statedim);             <span class="comment">% set initial covariance for PE proces noise</span>
0160 
0161 pNoise_PE.adaptMethod = <span class="string">'anneal'</span>;                     <span class="comment">% setup PE process noise adaptation method</span>
0162 pNoise_PE.adaptParams = [0.995 1e-7];                 <span class="comment">% We use the annealing method with a anneal factor of</span>
0163                                                       <span class="comment">% 0.98 and a variance floor of 1e-7</span>
0164 
0165 <span class="keyword">for</span> k=1:number_of_runs,
0166 
0167     fprintf(<span class="string">' [%d:%d] '</span>,k,number_of_runs);
0168 
0169     <span class="comment">% For dual estimation we iterate over the data, alternating between a state estimation step and a</span>
0170     <span class="comment">% parameter estimation step</span>
0171 
0172     <span class="keyword">for</span> j=2:N,
0173 
0174       <span class="comment">%--- First, we set the model parmaters of the state estimator using the surrent output of the parameter</span>
0175       <span class="comment">%--- estimator</span>
0176       InfDS_SE.model = InfDS_SE.model.setparams( InfDS_SE.model, Xh_PE(:,j-1), paramParamIdxVec);
0177 
0178       <span class="comment">%--- Now call the state estimator</span>
0179       [Xh_SE(:,j), Sx_SE] = <a href="../../.././ReBEL-0.2.7/core/srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>(Xh_SE(:,j-1), Sx_SE, pNoise_SE, oNoise_SE, noisy(:,j), [], [], InfDS_SE);
0180 
0181       <span class="comment">%--- And then the parameter estimator</span>
0182       [Xh_PE(:,j), Sx_PE, pNoise_PE] = <a href="../../.././ReBEL-0.2.7/core/srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>(Xh_PE(:,j-1), Sx_PE, pNoise_PE, oNoise_PE, noisy(:,j), [], Xh_SE(:,j-1), InfDS_PE);
0183 
0184     <span class="keyword">end</span>
0185 
0186     noisy_c = noisy(1:end-model.speech_taps+1);
0187     clean_c = clean(1:end-model.speech_taps+1);
0188     estim_1 = Xh_SE(1,1:end-model.speech_taps+1);
0189     estim_2 = Xh_SE(model.speech_taps,model.speech_taps:end);
0190 
0191     figure(1);clf; subplot(<span class="string">'position'</span>,[0.025 0.1 0.95 0.8]);
0192     p1 = plot(noisy_c,<span class="string">'k+'</span>); hold on;
0193     p2 = plot(clean_c,<span class="string">'b'</span>);
0194     p3 = plot(estim_1,<span class="string">'m'</span>);
0195     p4 = plot(estim_2,<span class="string">'r'</span>); hold off
0196     xlabel(<span class="string">'time'</span>);
0197     legend([p1 p2 p3 p4],<span class="string">'noisy'</span>,<span class="string">'clean'</span>,<span class="string">'estimate (0 lag)'</span>,<span class="string">'estimate (full lag)'</span>);
0198     title(<span class="string">'ReBEL Speech Enhancement Demo'</span>);
0199     axis tight
0200 
0201     figure(2);
0202     plot(Xh_PE'); hold off
0203     xlabel(<span class="string">'k'</span>);
0204     ylabel(<span class="string">'parameters'</span>);
0205     title(<span class="string">'Estimate of model parameters'</span>);
0206 
0207     drawnow
0208 
0209 
0210     mse(1,k) = mean((estim_1-clean_c).^2)/var(noisy_c);
0211     mse(2,k) = mean((estim_2-clean_c).^2)/var(noisy_c);
0212 
0213     fprintf(<span class="string">'  Normalized MSE : 0-lag estimate  = %4.3f   full-lag estimate =  %4.3f\n'</span>,mse(1,k),mse(2,k));
0214 
0215     <span class="keyword">if</span> dosound
0216       fprintf(<span class="string">'   Playing : noisy sample...'</span>);
0217       soundsc(noisy_c,8000,16);
0218       pause(1);
0219       fprintf(<span class="string">' 0-lag estimate...'</span>);
0220       soundsc(estim_1,8000,16);
0221       pause(1);
0222       fprintf(<span class="string">' full-lag estimate...'</span>);
0223       soundsc(estim_2,8000,16);
0224       pause(1);
0225       fprintf(<span class="string">' clean sample.\n'</span>);
0226       soundsc(clean_c,8000,16);
0227     <span class="keyword">end</span>
0228 
0229     <span class="comment">%-- Reset state estimates and covariance</span>
0230     Xh_SE(:,1) = zeros(InfDS_SE.statedim,1);
0231     Sx_SE = eye(InfDS_SE.statedim);
0232 
0233     <span class="comment">%-- Copy last estimate of model parameters to initial buffer position for next iteration...</span>
0234     Xh_PE(:,1) = Xh_PE(:,end);
0235 
0236 <span class="keyword">end</span>
0237 
0238 <span class="comment">%--- House keeping</span>
0239 
0240 <a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>(<span class="string">'../gssm'</span>);       <span class="comment">% remove relative search path to example GSSM files from MATLABPATH</span>
0241 <a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>(<span class="string">'../data'</span>);       <span class="comment">% remove relative search path to example data files from MATLABPATH</span></pre></div>
<hr><address>Generated on Tue 26-Sep-2006 10:36:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>